#This to obtain:
#Figure_2e.jpg
#Figure_2f.jpg


rm(list = ls())
library("sf")
library("ggplot2")
library("dplyr")
library("ggrepel")
library("raster")
library("RColorBrewer")
library("lemon")
library("sp")
library("ggsn")
library("gstat")
library("gtable")
library("lemon")

#setwd("C:/Users/bogdanp/")
setwd("/Users/bgpopescu/")

####################
#Reading Shapefiles#
####################

HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                layer="HRV_adm0_wgs")
HRV_adm0<-st_simplify(HRV_adm0, dTolerance = 0.005)

railroads_1879 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                          layer="railroads_1879")
railroads_1879$railroad<-"Railroads, 1879"
railroads_1879 = st_transform(railroads_1879, st_crs(HRV_adm0))
railroads_1879<-st_simplify(railroads_1879,  dTolerance = 0.005)

mil_col_short_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="krajna_line_GCS")
mil_col_short_polygon<-st_simplify(mil_col_short_polygon,  dTolerance = 0.005)


##############################
#Elevation for Croatia#
##############################

#Reading Rasters
euro_raster <- raster("./Dropbox/Legacies_Project/Analysis/data/euro_raster.tif")

HRV_adm0$y_lat<-sf::st_coordinates(st_centroid(HRV_adm0))[,1]
HRV_adm0$x_lon<-sf::st_coordinates(st_centroid(HRV_adm0))[,2]

max_lat<-max(HRV_adm0$y_lat) + 4
min_lat<-min(HRV_adm0$y_lat) - 4

max_lon<-max(HRV_adm0$x_lon) + 4
min_lon<-min(HRV_adm0$x_lon) - 4

#e <- extent(-160, 10, 30, 60)
e <- extent(min_lat, max_lat, min_lon, max_lon)
rc <- crop(euro_raster, e)	

#Lower resolution raster
rc.aggregate <- aggregate(rc, fact=1)
res(rc.aggregate)

#Get rid of really high mountains
rc <- reclassify(rc.aggregate, c(10000,Inf,NA))
rm(rc.aggregate)
rna <- reclassify(rc, cbind(NA, 0))
rm(rc)
#Convert rasters TO dataframes for plotting with ggplot
hdf <- rasterToPoints(rna); hdf <- data.frame(hdf)
colnames(hdf) <- c("X","Y","Hill")

#   Create vectors for colour breaks
b.hs <- seq(min(hdf$Hill),max(hdf$Hill),length.out=100)

map_road_graph<-ggplot()+
  geom_raster(data=hdf,aes(X,Y,alpha=Hill)) +
  scale_alpha(name = "Altitude in m", guide="none")  + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "black", linewidth = 0.8)+
  geom_sf(data = railroads_1879, colour = "red", linewidth = 0.8)+
  geom_sf(data = mil_col_short_polygon, fill = NA, colour = "blue", linewidth = 1.5)+
  scale_color_manual(name = "", values = c("Railroad, 1879" = "red"))+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(legend.position="left")+
  theme_bw()+
  #labs(x = "Longitude", y="Latitude")+
  #ggtitle("Distribution of Areas Covered by Zadrugas in 1895")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  ggsn::scalebar(x.min = 18.5, x.max = 19.1,
           #Above you mention how long the scalebar should be
           y.min = 46.4, y.max = 46.7,
           #Above you mention how thick the scalebar should be
           dist = 50, dist_unit = "km",
           transform = T, model = "WGS84",
           location = "topright",
           st.size = 4,
           #Above you have the font size of the numbers below the scalebar
           st.dist =0.2,
           #Above you have the distance between the bar and the text, as a proportion of the y axis.
           height=0.2)



ggsave(map_road_graph, 
       file="./Dropbox/Legacies_Project/Paper/figures/Figure_2f.jpg", 
       height=19.42, width=20.2,
       units = "cm", dpi=300)



#Reading Polygons
HABS <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                layer="all_hasburg_GCS")
sf::sf_use_s2(FALSE)
HABS<-st_simplify(HABS,  dTolerance = 0.005)


all_coasts<- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                     layer="only_coasts_GCS")
all_coasts<-st_simplify(all_coasts, dTolerance = 0.005)


map_road_graph2<-ggplot()+
  geom_sf(data = HABS, colour = "grey60")+
  geom_sf(data = railroads_1879, colour = "red", linewidth = 0.8)+
  geom_sf(data = mil_col_short_polygon, fill = NA, colour = "blue", linewidth = 1)+
  geom_sf(data = all_coasts, colour = "black", fill = NA, linewidth = 0.5)+
  scale_color_manual(name = "", values = c("Railroad, 1879" = "red"))+
  coord_sf(ylim=c(40.815009,53), xlim=c(8.556495, 27))+
  theme(legend.position="left")+
  theme_bw()+
  #labs(x = "Longitude", y="Latitude")+
  #ggtitle("Distribution of Areas Covered by Zadrugas in 1895")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  ggsn::scalebar(x.min = 10, x.max = 14,
                 #Above you mention how long the scalebar should be
                 y.min = 52, y.max = 53,
                 #Above you mention how thick the scalebar should be
                 dist = 200, dist_unit = "km",
                 transform = T, model = "WGS84",
                 location = "topright",
                 st.size = 4,
                 #Above you have the font size of the numbers below the scalebar
                 st.dist =0.2,
                 #Above you have the distance between the bar and the text, as a proportion of the y axis.
                 height=0.2)

ggsave(map_road_graph2, 
       file="./Dropbox/Legacies_Project/Paper/figures/Figure_2e.jpg", 
       height=19.42, width=20.2,
       units = "cm", dpi=300)
